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In many situations a BCS-type superconductor will develop an imbalance between the populations 
of the hole-like and electron-like spectral branches. This imbalance suppresses the gap. It has been 
noted by Gal'perin, Kozub and Spivak [6 that at large imbalance, when the gap is substantially 
suppressed, an instability develops. The analytic treatment of the system beyond the instability 
point is complicated by the fact that the Boltzmann approach breaks down. We study the short 
time behavior following the instability, in the coUisionless regime, using methods developed by 
Yuzbashyan et al. [T81IT9] . 



The excitation spectrum of a BCS superconductor consists of an electron-like and a hole-like branch. While the 
two are equally populated in equilibrium, non-equilibrium states may have a 'branch imbalance'. For example, if a 
superconductor is placed in an NSN junction as in the experiment of Clarke [U HI |3] , the injected quasi-particles can 
be primarily electron-like, and the electron-like branch more heavily populated than the hole-like branch in the steady 
state. Superconducting wires, where branch imbalance arises at phase-slip centres [4^, give another example. 

This imbalance suppresses the spectral gap [5]. If large enough, it returns the system to the normal state, which is 
known to be unstable below - this is the Cooper instability. It was realized by Gal'perin, Kozub and Spivak [SlIT] 
that the Cooper instability is a limiting case of a more general instability which afflicts any state beyond a critical 
level of imbalance and results in oscillations of the order parameter. Traditional approaches to the dynamics fail 
when the instability occurs, and for this reason the supercritical behavior of the superconductor is as yet unknown. 
It is important to resolve this, since situations in which large imbalance occurs are quite natural, for example in NSN 
junctions at large enough injection rate, or long superconducting wires held at sufficiently high voltage. 

This paper makes a step in this direction by treating the short-time dynamics in the limit where dissipative processes 
act slowly in comparison with the BCS dynamics (in particular the oscillation of the gap). We work in the regime 
Tc — T <ti Tc, when the slow relaxation of imbalance allows the system to reach a 'quasiequilibrium' whose deviation 
from equilibrium can be characterized only by the amount of imbalance. We also assume the initial conditions are 
close to the unstable stationary solution, in a sense defined below. 

Since the BCS dynamics are characterized by a time scale of order l/(Tc — T), the assumed separation of scales 
is ^/(Tc — T) <C Te, where is the time scale associated with energy relaxation. In a metallic superconductor with 
Debye energy 6 » and r, ~ Q'^/T^, this gives only the unrestrictive (T^ - T)/Tc > T^/Q^. 

Often one can avoid the complexities of the microscopic BCS dynamics of a superconductor with a simpler effective 
description such as the Ginzburg-Landau or the Boltzmann kinetic equation. Two time scales are important in 
deciding whether either is appropriate: the inelastic quasiparticle relaxation time, r^, and the time scale ta over 
which the order parameter varies significantly. 

When Te <C TA, the quasiparticle distribution rapidly reaches a local equilibrium characterized by the order pa- 
rameter A{r,t) — |A(r, <) |e*-^'^'''*\ with dynamics described by the Ginzburg-Landau equations for A(r,t). We are 
interested in the opposite limit, 3> ta, which is usually tackled with the Boltzmann kinetic equation [5] for the 
quasiparticle distribution function n{f,p,t): 



Here e is the energy of a quasiparticle state and the functional I{n} accounts for impurity scattering and collisions 
between electrons or between electrons and phonons. The kinetic equation must be supplemented with the self- 
consistency equation, which determines |A(f', i)| and thus the quasiparticle energies: 



I. INTRODUCTION 



dn de dn de dn 



I{n}. 



(1) 



dt dp df dr dp 




(2) 



Here S^j; = p^ /2m — /i, /i is the chemical potential, and A the BCS coupling constant. In addition there is a 'neutrality 
condition' involving the phase x of the order parameter. This condition arises from the continuity equation and ensures 
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the conservation of charge - it is independent of ([T]) in the case of a superconductor. In the spatially homogeneous 
case, if we define <I> = |^ + e<y9 (where Lp is the electric potential) and = + it takes the form 

"^"^^<P = <i>. (3) 



This integral quantifies the branch imbalance. The electron- and hole-like branches are distinguished by the sign of 
^. In equilibrium the two branches are identically populated since the quasiparticle energy (e^ = C'^ -I- |Ap) is even in 
this quantity, and $ vanishes. 

It is important that while branch imbalance is absent in equilibrium, its relaxation rate tq diverges as is 
approached [3] H]- ~ t^Tc/A. (This is in the absence of oscillations of the gap, and assuming relaxation due to 
electron-phonon collisions j3j[8].) Thus when <i> 7^ 0, the quasiparticle distribution reaches a 'quasiequilibrium' on a 
time of order Tj, characterized by distinct chemical potentials for the separately equilibrated hole-like and electron- like 
branches O O [HI HHj • This distribution is given below (10). Inserting it into the self-consistency equation ([t]) reveals 



that imbalance suppresses the gap relative to Aq, its value when <f> = 0: 

|Ap=A2-2$2. (4) 

At <& = Ao/-\/2, the gap is completely suppressed, and the system is returned to the unstable normal state. The 
instability appears earlier, at an intermediate value <I>c. 

The Boltzmann description cannot handle the system after the instability takes hold, as modes are excited in 
which Cooper pairs posses non-zero relative phases (the relative phases of the S-(C) = S2:(^) — *Sy(^), in the language 
of Anderson pseudospins Only the overall phase of the condensate is retained in the Boltzmann approach, 

effectively restricting the system to a subclass of solutions where the Cooper pairs precess in phase. 

In this situation we must return to the Gorkov equations describing the mean-field dynamics of the individual 
Cooper pairs T^. We study here only the limit in which dissipative processes are neglected (r^ 00) and the 
dynamics controlled purely by the BCS Hamiltonian. This was done for the Cooper instability by Barankov, Levitov 
and Spivak in [13], yielding a 'soliton train' of peaks in the gap. The problem was further discussed by Warner and 
Leggett in 15], by Barankov and Levitov in [THIII^, and by Yuzbashyan and Dzero in [IT. The integrability of the 
mean-field BCS system was established by Yuzbashyan, Altshuler, Kuznetsov and Enolskii in the papers [18, and 
a general framework for addressing its dynamics was developed, which we use here to tackle the more involved case 
of imbalance. We confine ourselves to the integrable BCS Hamiltonian because it can be treated analytically; it was 
however shown by Barankov and Levitov in ^D] for the case of the Cooper instability that the gap oscillations survive 
the breaking of integrability. 

Our results show the emergence of oscillatory behavior at the instability point. On much larger time scales (~ Tg) 
dissipation will modulate the form of these oscillations and determine the final fate of the system. Such an analysis 
is beyond the scope of this paper, but the solutions we present may be thought of as candidates to be found either 
at very large times after the onset of the instability or at intermediate times on the way to the asymptotic behavior, 
and are a first step in calculating the long-time behavior. We discuss this briefly in Section [V] 

The ability to 'switch on' the pairing interaction in ultracold trapped gases via tunable Feshbach resonances [HJ 
W2\ E51 [231 [25J [53] means that one might hope to observe the oscillations of the order parameter directly in the case of 
the Cooper instability [T31 [TH [121 [Sj dZ] . On the other hand, it is not clear whether one could create imbalance in a 
sufficiently controlled fashion to observe the instability it creates [33]. In contrast, imbalance can be created easily in 
metals, for example in tunnel junctions [I1[H[5]. Here the direct observation of the coUisionless dynamics is unlikely 
due to short dissipation times, and its significance is instead in its effect on processes at longer timescales such as the 
relaxation of imbalance. 

The structure of the paper is as follows. Section [ll] presents the equations of motion for the system and the relevant 
initial conditions, and includes the final result of the analysis. Sectio n |III| describes the formal solution to BCS 
dynamics [18j [THl [2H1 HSj , which is applied to the relevant case in Section IV Wc conclude in Section |vj 



II. GORKOV'S NONLINEAR EQUATIONS & SOLUTION WITH IMBALANCE 

The system of equations describing the BCS superconductor in the non-dissipative regime were derived by Volkov 
and Kogan [12] in the Keldysh Green's function formalism: 

siO = -?(Ox(2A„2A^,-2a, (5) 
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where s is defined by the following Keldysh Green's functions: 

.,(e) = ([cT(o,4(e)]) 



s-(e) = ([cT(0,ci(e)]> (6) 



iSy, and: 



A = A,-zAy= - / (7) 



This system of equations can be derived from a classical Hamiltonian: 

H = l2^s.{0di-j\A\'; {s^{0,s,{a}^e,,kSk{OS{^~a- (8) 

This is the mean-field BCS Hamiltonian written in terms of Anderson pseudo-spins, with an up (resp. down) spin 
representing a full (empty) Cooper pair. Singly-occupied pairs decouple from the order parameter dynamics (as can 
be seen from the BCS Hamiltonian which involves only pair-to-pair scattering) and correspond to zero-length spins. 
One may derive the Gorkov equations heuristically as a mean-field approximation to the BCS Hamiltonian. 

The current paper aims to present solutions of ^ for imbalanced superconductors, using the approach of |18j . 
Before doing so, let us recall the instability of the stationary solutions of ^ in the presence of imbalance which was 
pointed out in fB|. Generally, stationary solutions of (|5| have the form: 

= , ^(l-2n(g)) (9) 

= ^7^^^^(l-2n(e)) 

where n(^) is the quasi-particle distribution function, which in the presence of imbalance is|34| [31 [5] 

n{0 = . , ^ : X (10) 

/ ^(c+$)2+|Ap-j.sign(g+<i.) \ ^ ^ 



exp 



T 



which is a Fermi-Dirac distribution for the quasi-particles, with different chemical potentials for the hole-like and 
electron-like spectral branches implemented by the term <I>sign(^ -I- $) in the exponent. $ parametrizes the amount 
of imbalance in the system. The expression for n(5) is valid when |^ -|- $| ^ |A|. The self consistency condition ([t]) 
is satisfied if [5]: 

|Ap = A2-2*2 (11) 

where Aq is the order parameter at the temperature T in the absence of imbalance (for the case $ = 0). Thus 
imbalance between electron- like and hole-like excitations suppresses the order parameter I35|. 



One can easily see that the distribution is unstable when $ = Ao/-\/2, and A = by (11). At this point it becomes 



n(0 = ^ (12) 

exp(MS±5))+l 

which is nothing but the quasi-particle distribution of a normal Fermi gas in the excitation representation, where an 



artifical distinction between hole-like and electron-like excitations is made at <^ = — <i>. So, the peculiar form of (12 1 
is an artifact of the excitation representation, and it just describes a normal metal placed at T < Tc- The Cooper 
instability of this metal presents itself as a instability of the stationary solution of equation ([s]). 

At (f> = however the solution is stable, and represents the equilibrium superconducting state. There must therefore 
be an onset of instability at some finite $ intermediate between and Ao/V2. In order to find this point and in order 
to give a quantitative characterization of the instability, a linear stability analysis was performed around this solution 
in looking for the presence of an unstable mode e~"^*, Im(ti;) > 0. This leads to an integral equation for uj: 



(13) 
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FIG. 1: The soliton for the Cooper instabihty of a normal metal, with 7 = 0.5, and a pulse of oscillations of |A(t)| for <I> = 7, 
7 — 0.5, A = 2, showing the qualitative nature of the solution in the two regimes A = and 27 < A. The solution interpolates 
smoothly between these cases. 



Here ^ = ^ + <&. The solutions of this equation were found to be: 

ReM=2ci>, llmHI- ^((T-r,)-|A|). (14) 

TT 

The instabihty arises when the right hand side of this expression for |Im(tt')| becomes positive at ^c- When the 



imbalance completely suppresses the gap (14) agrees well with the usual Cooper instability case. In fact all we shall 



need in the following are the orders of magnitude of the the following quantities: 

A - s'^T,, 7 = ilm(a>) - s^T„ $ - sT, (15) 

where s = \J (Tc — T) /T^, and we have also defined the parameter 7, the rate of instability, an important parameter 
which will appear frequently below. The order of magnitude of the different quantities could just as well have been 
taken from the Cooper instability case, as they remain the same when the imbalance completely suppresses the gap. 
Namely these orders of magnitude are not sensitive to the exact form of the distribution function, but rather to gap 



suppression and may be viewed as consequences of Eq. ( 11 ) 



The main result of the following analysis is the time dependence of the order parameter following a perturbation 
from the unstable stationary solution. This is a train of 'solitons' each of the form: 

Ait) = A (^1 + Je-2'*(*-x)sech(27i)^ (16) 

where A (as opposed to A(t)) denotes the value of the order parameter, which we can take to be real and positive, in 
the stationary but unstable initial state. The solitons are separated in time by a period that increases logarithmically 
with the smallness of the perturbation, 

io = 7"'log(^) (17) 

in the notation used below, x is a real constant (see below) which does not affect the qualitative nature of the soliton. 

The behaviour of |A(f)| consists in fast (at frequency 2$) ocillations within an envelope given by 
|1 ± (27/A)sech(27f)|. When 27/ A ^ 1 (recall that as A we return to the finite-temperature normal metal) 
this formula gives a soliton of the same shape as that found in [13 for the order parameter oscillations associated 
with the Cooper instability. In the opposite limit where the instability is small and 27/ A <C 1, we find instead 
pulses of oscillations of |A(i)J during which the value of |A(t)| averaged over over a period of the fast oscillation rises: 
|A(t)U, ~ A(l + 7VA2sech^(27t)). 



III. ELEMENTS OF THE INTEGRABLE STRUCTURE 



The instability discovered by Gal'perin et al. shows that at large times dissipation may take the system to a new 
state which is very different to the initial stationary solution. A similar situation exists when a normal metal is placed 
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at T < Tc. For comparison we briefly review here the known results on this instability, the current paper expounding 
on these results to include the case where imbalance is present. 

A normal metal at T < is unstable against the development of a gap (this is the celebrated Cooper instability). 
The short term dynamics of the system, just after the instability takes hold, consist in an oscillatory behavior of 
A. These oscillations damp out at large times because of dissipation (this process is not described by the pure 
BCS Hamiltonian) and the material is left in the superconducting state. The oscillatory behavior of A consists of 
a 'soliton train' which can be found by solving ^ for a system placed near the normal metal state - namely with 
initial conditions (|9]), where Hp is the Fermi distribution function plus a small perturbation. This is discussed in 

[nnHHiiiB]. 

The soliton train describes the short-time (<g; t^) behavior of the order parameter following the onset of the Cooper 
instability. A similar analysis will be presented here for the case when the initial state is given by (9| and (10), i.e. 
a superconductor with branch imbalance near Tc- We study the behavior of the system after a small perturbation is 
added to n(^). This describes the short time behavior after the instability has taken hold. The approach also yields 
solutions where the perturbation is rather larger. We believe that these solutions may provide further intuition about 
the possible routes the system may take once dissipation effects are taken into account. 

In order to find these solutions we apply the formalism developed in [18i il9j for the dynamics of the mean field 
BCS system, which draws heavily on the theory of integrable systems. We shall establish the notation and present 
the main concepts of the derivation, referring the reader interested in further details to the original papers. In these 
papers the spectrum is treated as discrete - here we assume a continuous spectrum. 

An important object in the integrable structure is the Lax matrix of the system, a 2 x 2 matrix depending on a 
complex parameter u, given by: 

««)^f„^ff. fir,!>) (18) 



where: 



B*{u) -A{u) 
J 

For any u, the matrix C(u) is time dependent via s, but its eigenvalues are not - a key to integrability. These 
eigenvalues constitute an infinite set of constants of motion[36j. They are labeled v{u), and given by: 

v{u) = ±V- det Ciu). (20) 

The analytic structure of v{u) is as follows: in general v(u) has branch cuts, with square root behavior around the 
branch points, parallel to the imaginary u axis, as well as a jump discontinuity on the real axis, on the support of 
the spectrum. The branch points Ei satisfy v{Ei) = 0. Other important quantities are the zeros of Biu), dubbed 
B{ui) = 0. 

An important simplification takes place if one is interested only in the long time behavior of the system (but still 
at times <C r^). After an initial transient the system exhibits periodic or quasi-periodic behavior whose frequency 
is dictated by the branch points of v{u). The jump discontinuity is only relevant for the initial transient. The 
oscillatory behavior following the transient is captured by a simpler system containing only a finite number of spins 
s*-*-*, 1=1, ...,5 -t- 1, where g -I- 1 is the number of branch cuts in v{u). The system with a finite number of spins is 
integrable, with a similar integrable structure to the infinite system, integrals being replaced by sums; namely if 

2 s^') 



B{u) = J2 



9+1 



are substituted into (18) then the eigenvalues of the matrix are constants of motion. It is convenient to define 



P{u) = Yiii''^ ~ Ci)i ^iid ^Iso a degree 2(/ -I- 2 polynomial Q(u) with zeros at Ei, through v(u) = | \/ Q(u) /P{u). The 

appearance of Q{u) signals the relevance of the algebraic Riemann surface defined by the curve y{u) = \J Q(u) to 
this problem. To find the configuration of the finite number of spins at any given time one must know Q{u), which 
is independent of time, and the time-dependent quantities Ui(t), defined by B{ui{t)) = 0. 



FIG. 2: The three-sheeted Riemann surface with branch points and cycles labeUed. 



To find the dependence of the Uj on time, it is best to make use of the connection of integrable systems and 
Riemann surfaces. A central theme in the study of Riemann surfaces are the cycles, which are the non-trivial closed 
curves on the Riemann surface (those that cannot be smoothly shrunk to a point). They will be denoted by bk and 
hk, k = 1, ...,(?. These cycles are depicted in Fig. [2j Another mainstay of the theory of Riemann surfaces are the 
Abelian (meromorphic) differentials on the surface. The so-called differentials of the first kind, which are everywhere 
holomorphic, form a g dimensional vector space for which a basis is: 

= (22) 

The differentials may be integrated around the non-trivial cycles to assist in the following definitions: 

Wfc.i = ^ / ill, i^k.i^if T = uj^'^uj'. (23) 

A familiar construction for genus one Riemann surfaces, which have the topology of a torus, allows us to represent 
the surface as a rectangle with opposite edges identified. The rectangle is characterized by its aspect ratio, which is 
an invariant of the Riemann surface as well. In the case of genus one the aspect ratio of the rectangle turns out to 



be equal to —it from (23 1. When the genus is higher than one we encounter a matrix, r, which is a generalization of 
the number r of the genus-1 case. The rectangle with opposite sides identified also has an analogue for higher genera: 
it is replaced by 2g (real) dimensional volume in given by /{Z^ui + Z^uj'). This 2g dimensional volume is an 
analog of the 2 dimensional rectangle in the genus-1 case, in that there exists an invertible map taking sets of points 
on the Riemann surface into it. This is given by J{{u}) : /{Z^cu + Z^uj'), where is the space of sets of g 

points {ui, . . . ,Ug} (these are not ordered sets, so permutations are considered equivalent): 

J,{{ui\Ui)^j2 r '^^r (24) 
,:=i •'Po 

The space C^/(Z^(jj + Z^uj') is called the Jacobian. The contour of integration from the arbitrary initial point Pq 
to the point u can wind around any of the cycles any number of times, so as a mapping to it is only defined 
up to the addition of an element of the lattice Z^lo + Z^lo' . This is however enough to give a well-defined map to 
'U^ jiZpLij + Z^oj'). One can check that in the case of genus one the mapping takes a point on the Riemann surfaces 
and maps it onto a rectangle whose aspect ratio is —ir, the rectangle being represented by C/{Zuj + Zoj'). 

The concept of the Jacobian is particularly important in the solution of the problem because one can show that 
the zeros of B(u), which are denoted by Ui, satisfy the equation: 

J{{n^{t)}Ul) = (ci, C2, . . . ,cg -f 2it) (25) 
where we have now written the time dependence of Ui explicitly, while the Ci are defined by: 

9 „uj{t=0} 



E 

.7 = 1 



dui. (26) 

E2j 
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The roots Ei of Q{u) are listed for our case in (32 1. Considered as constants of integration for the dynamics of the g 
roots Ui(t) of B(u), the q are g free complex variables determining the initial values Ui(0). However, not all initial 
values are permissible, i.e. an arbitrary c will not correspond to a configuration of the spins; there are g constraints 
on c. 

The Ui together with the spectral curve Q{u) contain all the information needed to find the configuration of the 
spins. The problem of finding the Ui is thus the problem of inverting the map J. This is a solved mathematical problem 
with a long history, which goes by the name 'the Jacobi Inversion problem' [50]. The solution can be obtained in 
terms of the Riemann 0-function. We are interested here in the order parameter, whose logarithmic derivative in time 
is given by Ui . The explicit solution of the inversion problem for BCS dynamics is presented in P^E] , and gives the 
time dependence of the order parameter as: 



m = ^ E = C'exp (2di^{Lo'^)~^x - t/3t) -) ) ^ (27) 



((2a;^)-i(x-d)|r) 



provided that rj and d are given by: 



The frequency /3 can be written in terms of the roots of Q{u) as /3 = 2^,- Ei. 

IV. SOLUTION WHERE IMBALANCE IS PRESENT 



In the stationary (but unstable) imbalanced state, Sz{S,) and S-(0 sjce given by the expressions (|9| and (10 1. In 
this case, the eigenvalue v{u) is given by: 

v^{u-^)^{u^ + A^)G{uf. (29) 

Throughout we use A (as opposed to A{t)) to denote the value of the gap before the perturbation. G{u), which 



was defined in (14 1 because it appeared in the linear stability analysis, has appeared again as a common factor in 
det C{u) — —A^iu)— B{u)B* (u) . This is not a coincidence - the connection between roots of v{u) and modes present 
in the solution is explained further in |28j . Since v{u)^ has six roots we are led to consider the three spin problem, 
which according to the arguments of the previous section represents the dynamics of the order parameter after an 
initial transient. The polynomial Q{u) for the three-spin problem has the same roots as w(it)^, given thus by: 

Q{u) = {u^ + A^){u - $ - ij)^iu - $ + i"/)^. (30) 

(Here we have shifted u by 4> - from the equations of motion in [18] we see that this can be compensated by giving 
the order parameter an additional phase factor.) Q{u) has a very particular structure, related to the fact that it is a 
stationary solution. It has single roots only at ±iA and the rest of its roots are double roots. This insures that the 
Riemann surface given by Q{u) = y'^{u) is of genus 0. 

We now wish to perturb the initial conditions. Unless we fine-tune the perturbation to avoid doing so, we will lift 
the degeneracy associated with the double roots: the polynomial Q{u) will now have 6 single roots (which must still 
occur in complex conjugate pairs, since Q has real coefficients): 

Q{u) = {u^ + A2) {{u-<^>- ijf + S^) ((?/-$ -f ij)^ + 6^) . (31) 

For the integrals in (26 1, (28 1 we need the definitions (Fig. |2|: 



(£^0, Ei,E2,E3, Ei, E5) = (-iA, iA, <i> - + id,<i> + - iS,^ + + id, ^-ij-iS). (32) 

We neglect here (for example) the small shift in the position of the pair of roots around the origin: whereas the 
splitting of the double roots has a qualitative effect, such shifts have negligible effect on the solution, involving only 
small shifts in the parameters A, 7 and <&, and a small change in the overall rate at which the phase of the order 
parameter rotates. Also, while in general S can be complex, in the regime of interest to us the solution is not sensitive 
to the phase of S except through the value of c. In the following we treat 6 as real unless otherwise stated. 
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After the perturbation the Riemann surface is of genus two and A{t) is given by (27) with genus two hypereUiptic 
^-functions. The expression is quite formidable, yet certain features can be clarified without much analysis. Most 
notably, the solution is quasi-periodic, with quasi-periods which can be deduced straightforwardly from the general 
periodicity properties of 9 -functions together with the particular form the matrices lu and r take in this case. 

Before continuing to the analysis of the small perturbation case, wc note that if the perturbation is large enough it 
can lead to the appearance of new roots for Q{u), and to higher spin solutions (with higher genera) - this case is too 
general for us to say much about. 



We assume that the parameter S is small, as discussed in Section 
expressions for the matrices r and {2Lu'^)~^{x±d), which figure in i27\. We also expand in s = \/(Tc — T)/Tc 
the lowest order terms for each element. Then we have, to leading order (in practice wc must make sure lower order 
terms do not contribute): 



V| We then take the leading order in S of the 

taking 



^log-M^) 
*7log-M^) 





2* 



itl log M- 




27r ^"8 1 75A2 



ilog-^^) 



i + ¥log-^(^) 



{2.^r\x±d) = ( ± ^ log (t) , {^t ± ^) log-i (A) ) + (2.^) 



(33) 



We will return to the vector c, which requires further analysis. The period matrix r is seen to have very large and 
very small elements on the diagonal, diverging or vanishing with \og^^{5). Because of this the 6'- function is well 
approximated by trigonometric functions. First we use the modular invariance of 0-functions to trade in our period 
matrix for one whose elements are all of order log((5). This is done via the identity: 



lA —K 
—K ih 



-hy2 



which results in a 'transformed' r-matrix, 



^/h 



yi - tV^ 



i{A + K'^/h) in/h 
in/h i/h 



(34) 



Ttr = 



log 



4$^ 
<5A 

' 47 



log (^) 



log(^) 2 log 



27 
7r(f> 



1 1 
1 



(35) 



Once the ^-function has been written in terms of a r matrix with only large elements, its degeneration into trigono- 
metric functions is easily obtained from the definition of the ^-function in terms of an infinite sum. 



0{u\r)^ E 



(36) 



The dominant terms in the sum will come from the to close to the the stationary point of the real part of the exponent: 



Too = — (Imr) "'^Imu 



(37) 



(too is not necessarily a vector of integers). When to = (to, n) deviates much from toq the exponential becomes rapidly 
smaller because of the largeness of t, such that the sum is dominated by only a few exponential terms - a (hyper-) 
trigonometric polynomial. Because r is logarithmic in 5 and appears linearly in the exponentials, the sub-dominant 
terms in the sum will be suppressed by powers of 5. 



A. Recovering the Cooper Instability 

By setting $ = Ao/-\/2 we completely suppress the gap, returning the system to the normal-metal state. In this 
limit we should see the absolute value of the order parameter execute the train of cosh~^ solitons found in previous 
work on the Cooper instability |13j . This behaviour corresponds to a simpler two-spin solution. This degeneration 
into a system described by fewer spins (or lower-genus ^-functions) is expected whenever we close a branch cut on the 
Riemann surface for Q{u) (here that joining the roots at ±iA), so long as the initial conditions are such that there 
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is a Ui pinned at the resulting double root (here ui = 0) [TS]. This is the case when the perturbation is such that 
A <C i5, as is appropriate if the perturbation is most significant near the Fermi surface [37] . Note that here A 7^ A(0), 
since the former does not include the effect of the perturbation. With these initial conditions, an appropriate vector 
of constants c is given by 



27* 



log(^)-2i^(log(^)-2), i). (38) 



Once we have expressions for the thcta functions in terms of the 'large' r-matrix, we extract the asymptotes in the 
manner described above. Consider one 6*- function, say that in the denominator of (p7|. To begin with we find the 



stationary point toq = ("ic'^o) via (37 1, and see that while the value of uq changes with time, mo ~ —1 for small A. 
If we take only this m, we have a genus-one ^-function as expected on general grounds. The denominator degenerates 
similarly, and the ratio has a quasiperiod 

io = 7"'log(^y) (39) 
corresponding to the quasiperiod of the genus two 0- functions. Numerator and denominator each have a stationary 



value of n, which can differ from that given by (37 1 



no{m) = -{lm{T22)) ^lm{u2 + mTi2) (40) 

and the argument of each sum is a Gaussian in n whose width is fixed by the (2, 2) component of the transformed t 
matrix. 

For generic t we can ignore all but one n for the denominator, giving a single exponential, but when the stationary 
value of n is close to halfway between two integers two values of n are of comparable importance, yielding an expression 
for 9 in terms of a trigonometric function. (The numerator behaves similarly half a period later, but at these points 
the ratio is insignificantly small.) The end result is a train of solitons separated by to, each of the form (to leading 
order) : 

Mt) = — — (41) 
^ ' cosh {2-ft) ^ ' 

The vector c and the overall normalization are fixed in the foll owi ng way. The values of the ut (including ui = 0) 

A 



determine B{u) up to A(0): B{u) — ^^^u{u — ^2)- Equation (|30| gives Q{u) in the A ^ hmit. Writing Q{u) in 



terms of A{u) and B{u), we see that A{u) = j^u{u — wi){u ~ W2) for some Wi that are either both real, or conjugate to 



each other. Matching the coefficients of Q{u) written in terms of A{u) and B{u) with (30) gives a family of acceptable 
solutions for U2(0) and A(0), corresponding to different stages in the time evolution. Choosing a particular U2(0) 
allows us to integrate to get the c above, and also fixes A(0) and thus the normalization of our solution. In the above 
formula we have omitted an overall phase rotation which corresponds to a redefinition of the chemical potential. 

The form of the above soliton conforms exactly with previous results [T31 [TB] for the oscillations of the order 
parameter following a sudden turn-on of the BCS interaction. Interestingly, it also conforms exactly with the result 
of the next section, where we assume (5 ^ A, in the limit that A ^ 7. This is despite the difference in the relative 
sizes of A and 5 in the two cases, which implies very different initial conditions for the Ui. 

B. Main Case 

We now consider the case (5 <C 7, A. 

1. Initial Conditions for Small Perturbations 



An important difference between this case and the previous is the value of c, which ( 26 1 gives in terms of the initial 
positions of our variables Ui(t = 0). These are the zeroes of B{u) at time t = Q, which in the unperturbed case 
coincide with the zeroes of G{u) as one can show using the self-consistency equation. For a small perturbation, the 
Ui remain close to the branch cuts at <I> ± 17; let us call them u±. It is useful to define z± as the distance of u± from 
the centers of the branch cuts in units of 5, which may be complex: 



u+(t = 0) = <& -I- Z7 -I- u_(i = 0) = (f> - n -f z-5* 



(42) 
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The integrals ( 26 1 defining c can then be given through z± using the function I± defined as 



2± 



dw 



ITT 



I± = - y 2 , = V - arcsinh(z±). (43) 

Ji Vw +1 ^ 

The expression for Ci is then given by [38]: 

up to corrections suppressed by 7^/<i>^ and zS/'^. 

The form of the solution is sensitive to the values of ci and C2. But as noted above, these are not independent 
parameters. The necessary constraints can be found in the following way. We first find expressions for z± in terms of 
the perturbations to A{u) and B{u). Let uj be the position of a root of B{u) before the perturbation is added, and 
Up its position afterward. Expanding B — Bq + SB about the initial position of the root uj tells us that, due to the 
perturbation, u travels a distance given by 

(up - ui) ~ -dB{ui)/B'f,{ui). (45) 

Similarly we can expand Q cx A{u)'^ + B{u)B*{u), taking into account the fact that it has a double root to begin 
with, and the fact that before the perturbation Aq and Bq are related by Bq{u) = ^^^Aq{u). This yields both the 
width of the branch cut (i.e. 26 or 26*) that is opened up by the perturbation, and the location of its centre, in terms 
of 6B{uj) and 6A{ui). We omit these formulas. Then the z± are given by 

(up — ui) — (displacement of centre of branch cut) ^^^^ 
(complex half-length of branch cut) 

Defining 6B{u) — ^6A{u) + 6C{u), all the 6 As and 6Bs disappear in favor of 6Cs, and we can expand without 
worrying about the relative size of 6 A versus 6B: 



^±»7 6C{^±ijl 

'^^^^6cw±wy ^ ^ 

Since 6C*{^ + ij) — 6C{^ — ij)* , this tells us that = <i)2/A^ to leading order, and that[3n] argz+/z_ is of 



order 7/$. Eq. (47 1 yields the following constraints on a or z±: 



Pi =Re (ci^^-cs*) = ^log4|z+z_| =-log^, (48) 

1 z_ 

P2 = Re (7C2) = - arg — ~ 0. (49) 

2 z+ 

These combinations of ci and C2 are precisely those necessary for the correct expansion of the ^-functions - for example 
Pi dictates which integers (m, n) give the leading order contributions to the representation of the theta function as a 



sum (36 1 



Having derived these results by expanding B{u), A{u) and Q{u), we must ask when they are valid. Assuming that 
6C{ui) 1 6C* {ui) is approximately of order one, we find that the roots of B{u) move a distance of order 6^/ A. This 

quantity must be much smaller than 7, the scale on which our initial polynomials vary. So a necessary condition for 

the validity of these approximations is 

(5<7A/$. (50) 
This excludes of course the Cooper instability case, where one root of B{u) is a distance of order $ from the start-points 



of the integrals in (43 1. Since B{u) vanishes as A —> 0, in this case it is not legitimate to assume that 6B <C Bq. 



2. Form of the Solution 



Again we use the 'transformed' 0-functions, and extract the dominant exponentials from the sums defining them, 
(36 1. The precise values of c depend on the nature of the perturbation, but the information obtained above is enough 
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to establish the nature of the solution up to (a) an overall shift in time and (b) a shift of the fast oscillations within 
their envelope. Up to such a shift, each soliton has the form (we neglect 7/<i> corrections): 

A{t) = A (^1 + ^e-2'**sech(27t)^ (51) 

and solitons occur at intervals of to — log(47/(5). 

More explicitly, taking into account the expressions for the shifts in terms of c, the first soliton has the form 

A{t) = A ^1 - ^exp ^-2i$t- C2* + ci$2 +log ^ + iarg5^sech^27(t - to/2) - hc2^^ . (52) 

From ( 43|47 1 , Im C2 is of order 7^ ^ log(<i> /A) , so that if 6 is sufficiently small the first soliton takes about half a period 
to appear. 



V. CONCLUSION 



We have found the short-time behavior of a BCS superconductor following a small perturbation to the imbalanced 



initial conditions given by Eq. 10 These initial conditions show a suppression of the gap |5] with increasing imbalance 
<&, and an instability ^ when $ > <i>c which becomes the celebrated Cooper instability when the gap is fully suppressed. 
As is increased beyond <I>c, the gap oscillations following upon the instability grow in magnitude. They take the form 
of a train of solitons, each of duration 7"^ and magnitude 7, and containing oscillations on the shorter timescale 
(Eq. 51 and pictured in Fig. [2]). These oscillations should be observable if appropriate initial conditions can be 



prepared in a controlled fashion. 

A stronger motivation for the work is that the oscillatory behavior is relevant to evolution on longer time scales 
{'^ Te) in experimental situations with large imbalance. In particular, the oscillations are relevant to the relaxation of 
the imbalance, which in the absence of the instability occurs at a rate which vanishes with the gap as T — > Tc [3j[9]. 
Understanding the short-time dynamics of (|5]) is a first step; to determine quantitatively what happens on long time 
scales it is necessary to compute how collisions modulate them. The moduli of the solution, i.e. the variables used 
to parameterize the kinds of short time behavior, can be taken to be the roots of Q{u). These roots or moduli vary 
slowly with time on account of collisions. The non-equilibrium state of the system at long times may correspond for 
example either to an unchanging set of moduli, or to a limit cycle in moduli space. Such an analysis is beyond the 
scope of the current paper, but will involve (51 ) and possibly generalizations. 

Our explicit expressions for the behavior of the order parameter apply when 5 <C 7. In this limit, where the 
solitons are widely spaced, the expressions simplify greatly; but the generalization to larger 5 may be necessary to 
treat imbalance relaxation in realistic situations. As an idealized Gendankenexperiment, our limit can be realized 
by instantaneously injecting electrons at the Fermi level to a system at the instability point - it can be shown using 
the definition (20 1 that such a perturbation increases the instability rate 7 while hardly increasing 5|40]. (A similar 
analysis shows that in a system with tunable interaction strength, increasing the coupling of a system with ^ <i>c 
does the same.) The subsequent evolution of such systems on time scales may increase 5 further. 

The present work neglects spatial inhomogeneities. Whether they change the picture qualitatively in systems larger 
than the coherence length remains to be investigated. Gap oscillations can also parametrically excite inhomogeneous 
modes, as shown in [51] for the Cooper instability. 

We have already mentioned that thermal processes act by slowly perturbing the dynamics considered. We have not 
mentioned thermal fluctuations in the initial conditions, which were analysed for the Cooper instability in [l51 132] . 
These fluctuations disappear when the level spacing, or the effective level spacing in a coherence length, goes to zero, 
but will cause variations in the parameters of the solution (e.g. our 5) when it is finite. While the evolution is of 
the same form in each realization, it was shown that variations in the parameters between different realizations are 
qualitatively important for averages (e.g. of the absolute value of the gap) over realizations. Such averages would be 
relevant to experiments involving direct observation of gap fluctuations. 
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